{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "WARNING: user expression has not supplied value_shape method or an element. Assuming scalar element.\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "d3d31a03cbf94404b6d17c6f63eb0436",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Plot(antialias=3, axes=['x', 'y', 'z'], axes_helper=1.0, background_color=16777215, camera=[0.0708601181276387…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "\"\"\"Compute the magnetic field B in an iron cylinder, the copper wires, and the surrounding vacuum.\"\"\"\n",
    "# https://fenicsproject.org/pub/tutorial/html/._ftut1015.html\n",
    "from fenics import *\n",
    "from mshr import *\n",
    "from math import sin, cos, pi\n",
    "\n",
    "a = 1.0   # inner radius of iron cylinder\n",
    "b = 1.2   # outer radius of iron cylinder\n",
    "c_1 = 0.8 # radius for inner circle of copper wires\n",
    "c_2 = 1.4 # radius for outer circle of copper wires\n",
    "r = 0.1   # radius of copper wires\n",
    "R = 2.5   # radius of domain\n",
    "n = 5     # number of windings\n",
    "\n",
    "# Define geometry for background\n",
    "domain = Circle(Point(0, 0), R)\n",
    "\n",
    "# Define geometry for iron cylinder\n",
    "cylinder = Circle(Point(0, 0), b) - Circle(Point(0, 0), a)\n",
    "\n",
    "# Define geometry for wires (N = North (up), S = South (down))\n",
    "angles_N = [i*2*pi/n for i in range(n)]\n",
    "angles_S = [(i + 0.5)*2*pi/n for i in range(n)]\n",
    "wires_N = [Circle(Point(c_1*cos(v), c_1*sin(v)), r) for v in angles_N]\n",
    "wires_S = [Circle(Point(c_2*cos(v), c_2*sin(v)), r) for v in angles_S]\n",
    "\n",
    "# Set subdomain for iron cylinder\n",
    "domain.set_subdomain(1, cylinder)\n",
    "\n",
    "# Set subdomains for wires\n",
    "for (i, wire) in enumerate(wires_N):\n",
    "    domain.set_subdomain(2 + i, wire)\n",
    "for (i, wire) in enumerate(wires_S):\n",
    "    domain.set_subdomain(2 + n + i, wire)\n",
    "\n",
    "# Create mesh\n",
    "mesh = generate_mesh(domain, 64)\n",
    "\n",
    "# Define function space\n",
    "V = FunctionSpace(mesh, 'P', 1)\n",
    "\n",
    "# Define boundary condition\n",
    "bc = DirichletBC(V, Constant(0), 'on_boundary')\n",
    "\n",
    "# Define subdomain markers and integration measure\n",
    "markers = MeshFunction('size_t', mesh, 2, mesh.domains())\n",
    "dx = Measure('dx', domain=mesh, subdomain_data=markers)\n",
    "\n",
    "# Define current densities\n",
    "J_N = Constant(1.0)\n",
    "J_S = Constant(-1.0)\n",
    "\n",
    "            \n",
    "# Define magnetic permeability\n",
    "class Permeability(UserExpression):\n",
    "    def __init__(self, markers, **kwargs):\n",
    "        self.markers = markers\n",
    "        super().__init__(**kwargs)\n",
    "    def eval_cell(self, values, x, cell):\n",
    "        if self.markers[cell.index] == 0:\n",
    "            values[0] = 4*pi*1e-7 # vacuum\n",
    "        elif self.markers[cell.index] == 1:\n",
    "            values[0] = 1e-5      # iron (should really be 6.3e-3)\n",
    "        else:\n",
    "            values[0] = 1.26e-6   # copper\n",
    "mu = Permeability(markers, degree=1)\n",
    "\n",
    "# Define variational problem\n",
    "A_z = TrialFunction(V)\n",
    "v = TestFunction(V)\n",
    "a = (1 / mu)*dot(grad(A_z), grad(v))*dx\n",
    "L_N = sum(J_N*v*dx(i) for i in range(2, 2 + n))\n",
    "L_S = sum(J_S*v*dx(i) for i in range(2 + n, 2 + 2*n))\n",
    "L = L_N + L_S\n",
    "\n",
    "# Solve variational problem\n",
    "A_z = Function(V)\n",
    "solve(a == L, A_z, bc)\n",
    "\n",
    "# Plot solution\n",
    "from vtkplotter.dolfin import plot\n",
    "plot(A_z,\n",
    "     #isolines={'n':20, 'lw':1.5, 'c':'black'} #not yet working\n",
    "    )\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
